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Abstract. A method for the identification of small inhomogeneities 
from a surface data is presented in the framework of an inverse scat- 
tering problem for the Helmholtz equation. Using the assumptions of 
smallness of the scatterers one reduces this inverse problem to an iden- 
tification of the positions of the small scatterers. These positions are 
found by a global minimization search. Such a search is implemented 
by a novel Hybrid Stochastic-Deterministic Minimization method. The 
method combines random tries and a deterministic minimization. The 
effectiveness of this approach is illustrated by numerical experiments. 
In the modeling part our method is valid when the Born approximation 
fails. In the numerical part, an algorithm for the estimate of the number 
of the small scatterers is proposed. 



1 Introduction 

In many applications it is essential to find small inhomogeneities from surface 
data. For example, such a problem arises in ultrasound mammography, where small 
inhomogeneities are cancer cells. Current X-ray mammography will be replaced by 
the ultrasound one because X-ray mammography has a high probability of creating 
new cancer cells in a woman's breast in the course of taking the mammography test. 
Other examples include the problem of finding small holes and cracks in metals and 
other materials, or the mine detection. The scattering theory for small scatterers 



originated in the classical works of Lord Rayleigh. It was developed in [15(] and [16 



where analytical formulas for the scattering matrix were derived for the acoustic and 
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electromagnetic scattering problems. In [15| and [17] inverse scattering problems 



for small bodies are considered. In 18 the problem of identification of small sub- 



surface inhomogeneities from surface data was posed and its possible applications 



were discussed. In [11] the results of the numerical experiments are presented for 
the problem of finding one small inhomogeneity from surface data. In Q a method 
for finding small inhomogeneities from tomographic data is proposed. The main 
result of our paper is the new optimization procedure which can be used for actual 
finding the small subsurface inhomogeneities from surface scattering data. Our nu- 
merical results demonstrate that the method proposed in this paper is potentially 
of practical importance. 

Consider a point source y of monochromatic acoustic waves on the surface of 
the earth. Let u(x,y,k) be acoustic pressure at the point x, and k > be the 
wavenumber. The governing equation is: 

[V 2 + fc 2 + k 2 v{x)} u = -5{x - y) in R 3 , (1.1) 

u satisfies the radiation condition at infinity, and v(x) is the inhomogeneity in 
the velocity profile, x = {x\,X2 1 x 3 ). 

Let us assume that v(x) is a bounded function vanishing outside of the domain 
D = l4f =1 An which is the union of M small nonintersecting domains D m , all 
of them are located in the lower half-space R?_ = {x : x% < 0}. Smallness is 
understood in the sense kp <C 1, where p :— ^ maxixjn^/jdiam-Djn}, and diam 
D is the diameter of the domain D. Practically kp ^ 1 means that kp < 0.1, in 
some cases kp < 0.2 is sufficient for obtaining acceptable numerical results. The 
background velocity in (1.1) equals to 1, but we can consider the case of fairly 
general background velocity 
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Denote z m and v m the positions of the gravity center of D m and the total 
intensity of the m-th inhomogeneity v m := J D v{x)dx. Assume that v m ^ 0. 
The inverse problem to be solved is: 

IP: Given u(x,y,k) for all {x,y) on P := {x : X3 = 0} at a fixed k > ; find 
the number M of small inhomogeneities, the positions z m of the inhomogeneities, 
and their intensities v m . 

2 Method of solution 

Let us introduce the following notations: 

P := {x : x 3 = 0}, (2.1) 

{x 3 ,Vj}:=^, l<j<J, Xj , y, G P, (2.2) 
are the points at which the data u(xj,yj,k) are collected, 

fc > is fixed, (2.3) 



. , , exp(ifc|a; — y\) ,„ 

g(x,y,k :=— ^ f^, 2.4 

4tt\x — y\ 



Gj(z) :=G(£j,z) := g(x J ,z,k)g{y j ,z,k), 



(2.5) 
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To ■— k 2 ' 



M 2 



fj - Gj(z m )v r , 



(2.7) 



<S>(zi, . . . , Zm, vi,... ,vm) ■= ^2 

j—l m—l 

The proposed method for solving the (IP) consists in finding the global mini- 
mizcr of function (2.7). This minimizer {z\, . . . , zm, v±, . . . , %) gives the estimates 
of the positions z m of the small inhomogeneities and their intensities v m . 

The above approach can be justified as follows. The Helmholtz equation (1.1) 
with the radiation condition is equivalent to the integral equation 

M 

u(x 1 y,k)=g(x,y,k) + k 2 y2 g{x, z,k)v(z)u(z,y,k)dz. (2.8) 

m=1 JD m 

For small inhomogeneities the integral on the right-hand side of (2.8) can be 
approximately written as 

g(x,z,k)v(z)u(z,y,k)dz w / g(x, z,k)v(z)g(z,y,k)dz 

D m J D m 

~G{x,y,z m ) I vdz := G%z m )v m , 1 < m < M, (2.9) 

where z m is a point close to z m , and u(z, y, k) under the sign of the integral in 
(2.8) can be replaced by g(x,y, k) with a small error e 2 provided that 

e := c M(kp) 2 < 1, (2.10) 

where 

p= max p m , p m :=\diamD m , c := max \v(x)\, (2.11) 

l<m<M 2 £CGR 3 

and M is the number of inhomogeneities. 

Note that the sufficient condition for the validity of the Born approximation 
for equation (2.8), is the smallness of the norm in L 2 {D) of the integral operator 
in (2.8). The above condition can be written as 

Mcak 2 p 2 « 1. (*) 

If x, y G P, then 

\u-g\< 0(Mk 2 c a p 3 d- 2 ), 
where d > is the minimal distance from D m to P. 

The relative error of the first approximation in (2.9) for x, y <E P is 0(Mk 2 cop 3 d~ 1 ), 
where d > is the minimal distance from D m to P. If d > is not too small, 
cop 3 := V is not too large, then the above error is small if 

ei := Mk 2 c p 3 d~ 1 « 1, (**) 

where ei is a dimensionless quantity. Condition (*) can be written as Mk v « 
1. This condition is much stronger than condition (**), which can be written as 

<< 1, since d » p. Therefore the approximation (2.9) is applicable when the 
Born approximation may fail. 

If the condition 

c a p 3 = const as p — > (* * *) 
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holds, then condition (*) for the validity of the Born approximation is violated 
because Mk p v — ► oo as p — » 0. 

Condition (***) corresponds to the scattering by delta-type inhomogeneities. 
For such scattering the Born approximation is not applicable, but a modified rep- 
resentation of the type (2.9) is valid (see Q, p. 113, formula (1.1.33), and [g], p9| ). 

From (2.9) and (2.6) it follows that 

M 

h ~ X! Gj{z m )v m , Gj(z m ) := G(£j,z m , k) (2-12) 

m— 1 

Therefore, parameters z m and u m can be estimated by the least-squares method 
if one finds the global minimum of the function (2.7): 

®(zi, . . . , zm, v%,..., v M ) = min . (2-13) 

The function $ depends on M unknown points z m G Rl, and M unknown 
parameters v m , 1 < m < M. The parameter M, the number of the small inho- 
mogeneities, is estimated by the algorithm described in Sec. 3 below, see also a 
discussion in the beginning of Section 4. 

3 Hybrid Stochastic-Deterministic Method 

Let the inhomogeneities be located within a box 

B = {(xi,X2,xs) : — a < X\ < a, —b<x 2 <b, < x 3 < c} , (3.1) 
and their intensities satisfy 



<v m < v max . (3.2) 

Then, given the location of the points Zi,Z2,...,Zmi the minimum of $ in 
(2.13) with respect to the intensities v\, t>2, . . . , dm can be found by minimizing the 
quadratic function in (2.7) over the region satisfying (3.2). This can be done using 
normal equations for (2.7) and projecting the resulting point back onto the region 
defined by (3.2). Denote the result of this minimization by <t, that is 

$(zi, z 2 , . . . , z M ) =m.m{3>(z 1 ,z 2 ,...,ZM,Vi,V2,...,v M ) ■ ^ 
< v m < v max , 1 < m < M} 

Now the minimization problem (2.13) is reduced to the 3A/-dimensional re- 
strained minimization 

®(z!,z 2 , ■ • • ,z M ) = min, z m eB, 1 < m < M . (3.4) 

Note, that the dependency of $ on its 3M variables (the coordinates of the 
points z m ) is highly nonlinear. In particular, this dependency is complicated by 
the computation of the minimum in (3.3) and the consequent projection onto the 
admissible set B. Thus, an analytical computation of the gradient of $ is not 
computationally efficient. We have used Powell's quadratic minimization method 
to find local minima. This method uses a special procedure to numerically ap- 
proximate the gradient, and it can be shown to exhibit the same type of quadratic 
convergence as conjugate gradient type methods (see 

In addition, the exact number of the original inhomogeneities M or i g is unknown, 
and its estimate is a part of the inverse problem. In the HSD algorithm described 
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Figure 1 Objective function &(z r ,Z2, 23, 24, 25, 26) , — 2 < r < 2 

below this task is accomplished by taking an initial number M sufficiently large, so 
that 

M ong < M , (3.5) 

which, presumably, can be estimated from physical considerations. After all, our 
goal, is to find only the strongest inclusions, since the weak ones cannot be distin- 
guished from the background noise. The Reduction Procedure (see below) allows 
the algorithm to seek the minimum of $ in a lower dimensional subsets of the 
admissible set B, thus finding the estimated number of inclusions M. Still an- 
other difficulty in the minimization is a large number of local minima of $. While 
this phenomena is well known for the objective functions arising in various inverse 
problems, we illustrate this point on Figure 1. 

The original configuration (z\, £2, • • • , zm, vi,V2, - ■ ■ , %) , M orig = 6 is pre- 
sented in Table 1. Figure 1 shows the values of the function $(z r , z 2 , z 3 , £4, Z5, Ze), 
where 



z r = (r, 0, 0.520) , -2 < r < 2 



6 



S.Gutman and A.G.Ramm 



and 

z 2 = (-1,0.3,0.580). 

The plot shows multiple local minima and almost flat regions. 

A direct application of a gradient type method to such a function would result 
in finding a local minimum, which may or may not be the sought golobal one. In 
the example above, such a method would usually be trapped in a local minimum 
located at r — —2, r = —1.4, r = —0.6, r = 0.2 or r = 0.9. While the desired global 
minimum at r — 1.6 will be found only for a sufficiently close initial guess 1.4 < r < 
1.9. While various global minimization methods are known (see below), we found 
that the most efficient way to accomplish the task for this Inverse Problem was to 
design a new method (HSD) combining both the stochastic and the deterministic 
approach to the global minimization. Deterministic minimization algorithms with 
or without the gradient computation, such as t he conjugate gradient methods, are 
known to be highly efficient (see J|, |(| 10|). However, the initial guess should 



be chosen sufficiently close to the sought minimum. Also such algorithms, as we 
mentioned above, tend to be trapped at a local minimum, which is not necessarily 
close to a global one. A new deterministic method is proposed in [§] and 
which is quite efficient according to J3|. On the other hand, various stochastic 



are more 



minimization algorithms, e.g. the simulated annealing method [12, 13| 
likely to find a global minimum, but their convergence can be very slow. We 
have tried a variety of minimization algorithms to find an acceptable minimum of 
(3.4) or (2.13). Among them were the Levenberg-Marquardt Method, Conjugate 
Gradients, Downhill Simplex, and Simulated Annealing Method. None of them 
produced consistent satisfactory results. 

Among the minimization methods combining random and deterministic searches 



we mention Deep's method fa] and a variety of clustering methods [20], [21|. An 



application of these methods to the particle identification using light scattering is 



described in [22|. The clustering methods are quite robust but, usually, require a 
significant computational effort. The HSD method is a combination of a reduced 
sample random search method with certain ideas from Genetic Algorithms (see e.g. 
|)|). It is very efficient and seems especially well suited for low dimensional global 
minimization. Further research is envisioned to study its properties in more detail, 
and its applicability to other problems. 

While the steps of the Hybrid Stochastic-Deterministic (HSD) method are out- 
lined below, its basic idea is to start with a random search for the minimum. If at 
a certain configuration x s the minimized function is judged to be sufficiently small, 
this configuration is a candidate for the initial guess of a deterministic minimization 
method. However, this minimization is applied only after the points with smaller 
intensities in x s are dropped (Step 2 below), and all sufficiently close points are 
eliminated (Step 3). Now a deterministic minimization is performed in the result- 
ing lower dimensional subspace. The minimizer Xd is likely to represent at least 
some of the sought inhomogeneities. These points are supplemented with randomly 
chosen ones to obtain a full configuration (Step 1). In this manner the iterations 
of the random searches and the deterministic minimization continue till a tolerance 
criterion is satisfied. As the algorithm finds configurations with smaller and smaller 
values of the objective function, the likelihood of the finding the global minimum 
increases, while the irrelevant local minima are eliminated. 

Let us outline the steps of a possible implementation of the HSD algorithm 
for the surface data inverse scattering problem described in Sections 1 and 2. We 
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assume that all the inhomogeneities are located in a (known) box B, and M is 
chosen according to (3.5). Our goal is to locate the points ,5n, N < M 

which minimize $ in (3.4). 

Hybrid Stochastic-Deterministic (HSD) Method. Let Po, T max , e s , ej, e^, 

and e be positive numbers. Let N = 0. 

1. Generate randomly M — N additional points zy+i, . . . , zm to obtain a full 
configuration z\, . . . , zm- Compute P s = $(zi, %2 ■ ■ ■> zm)- Save the result- 
ing best fit intensities v m , 1 < m < M (see 3.3). If P s < Pq€ s then go to 
step 2, otherwise repeat this step 1. 

2. Discharge all the points with the intensities v m satisfying v m < v max Ei. Now 
only N < M points z±, z<i . . . , Zjy remain in the box B. 

3. If any two points z m , z n in the above configuration satisfy \z m — z n \ < td,D, 
where D = diam(B) 1 then eliminate point z n , change the intensity of point 
z m to v m + v ni and assign N := N — 1. This step is repeated until no further 
reduction in N is possible. 

4. Run a restrained deterministic minimization of $ in 3N variables, with the 
initial guess at the configuration determined in step 3. Let the minimize! - 
be si\, ... , zn- If P = ^K^i) • ■ ■ i zn) < c then save this configuration, and go 
to step 6, otherwise let Po = P, an d proceed to the next step 5. 

5. Keep intact N points Z\,... ,zn- If the number of random configurations 
has exceeded T max (the maximum number of random tries), then save the 
configuration and go to step 6, else go to step 1, and use these N points 
there. 

6. Repeat steps 1 through 5 n max times. 

7. Find the configuration among the above n max ones, which gives the smallest 
value to $. This is the best fit. 

Step 1 is the stochastic part of the algorithm. In step 2 the points with low 
intensities are discharged, since they, most likely, are artifacts contributing on the 
level comparable with the background noise. Step 3 is the Reduction Procedure 
replacing two nearby inclusions with one of the joint intensity. Steps 2 and 3 
lower the dimensionality of the minimization domain, thus greatly reducing the 
computational time needed to perform the deterministic minimization of the step 
4. We have used Powell's minimization method (see jl for a detailed description) 
for the deterministic part, since this method does not need gradient computations 
and converges quadratically near quadratically shaped minima. Also, in step 1, an 
idea from the Genetic Algorithm's approach |^] is implemented by keeping only the 
strongest representatives of the population and allowing a mutation for the rest. 

4 Numerical results 

The algorithm was tested on a variety of configurations. Here we present the 
results of just two typical numerical experiments illustrating the performance of 
the method. The data was simulated according to (2.12). In both experiments the 
box B is taken to be 

B = {(xi,x 2 ,x 3 ) : —a < xi < a, — b < x 2 < b, < x 3 < c} , 

with a = 2, 6=1, c=l. The frequency k = 5, and the effective intensities v m are 
in the range from to 2, (see Section 1). The values of the parameters defined in 
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Table 1 Actual inclusions. 



Inclusions 




X2 


^3 


V 


1 


1.640 


-0.510 


0.520 


1.200 


2 


-1.430 


-0.500 


0.580 


0.500 


3 


1.220 


0.570 


0.370 


0.700 


4 


1.410 


0.230 


0.740 


0.610 


5 


-0.220 


0.470 


0.270 


0.700 


6 


-1.410 


0.230 


0.174 


0.600 



Section 3 were chosen as follows 

Po = 1 , T max = 1000, e s = 0.5, e, = 0.25, e d = 0.1 , e = 10" 5 , n max = 6 

In both cases we searched for the same 6 inhomogeneities with the coordinates 
xi,X2,xs and the intensities v shown in Table 1. 

Parameter M was set to 16, thus the only information on the number of in- 
homogeneities given to the algorithm was that their number does not exceed 16. 
This number was chosen to keep the computational time within reasonable limits. 
Still another consideration for the number M is the aim of the algorithm to find 
the presence of the most influential inclusions, rather then all inclusions, which is 
usually impossible in the presence of noise and with the limited amount of data. 

Experiment 1. In this case we used 12 sources and 21 detectors, all on the 
surface X3 = 0. The sources were positioned at {(—2 + 0.333 + 0.667i, — 0.5 + 
l.Oj, 0), i = 0, 1, . . . , 5, j = 0, 1}, that is 6 each along two lines X2 = —0.5 and 
xi = 0.5. The detectors were positioned at {(—2 + 0.667«, — 1.0 + 1.0,7,0), i — 
0, 1, . . . , 6, j = 0, 1, 2}, that is seven detectors along each of the three lines X2 = 
— 1, X2 = and X2 = 1. This corresponds to a mammography search, where the 
detectors and the sources are placed above the search area. The results of the 
identification are shown in Tables 2 and 3 for different noise levels S in the data. 
See Figure 2. Table 3 has only 4 lines showing that the program has identified the 
presence of only 4 indicated inclusions, while missimg the other 2. 

To evaluate the performance of the algorithm, we have made 10 independent 
runs of the program for each of the noise levels S = 0.00, 5 = 0.02, and S = 0.05. 
In the noiseless case all 6 original inhomogeneities were identified every time. The 
program performed 2 or 3 deterministic minimizations in each run, spending from 
30 to 80 percent of time in this minimization. For 5 — 0.02 the perfect identification 
also was obtained in every run. The deterministic minimization was performed 2-4 
times, spending there 10 to 25% of the computational time. For 5 = 0.05 out of 
10 runs the program found 4 and 5 inhomogeneities in one run each, and all 6 
inhomogeneities in the rest 8 runs. The deterministic minimization was performed 
from 1 to 3 times in each run for the total of 5 to 20% of the computational time. 

Experiment 2. In this case we used 8 sources and 22 detectors, all on 
the surface x% = 0. The sources were positioned at {(—1.75 + 0.5i, 1.5,0), i = 
0,1,..., 7,. 7 = 0,1}, that is all 8 along the line x 2 = 1.5. The detectors were 
positioned at {(-2 + OAi, 1.0 + 1.0,7,0), i = 0, 1, . . . , 10, j = 0, 1}, that is eleven 
detectors along each of the two lines X2 = 1 and X2 = 2. This corresponds to a mine 
search, where the detectors and the sources must be placed outside of the searched 
ground. The results of the identification are shown in Tables 4 and 5 for different 
noise levels S in the data. See Figure 3. Because the sources and the detectors 
are further away from the inhomogeneities than in Experiment 1, this presents a 
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Table 2 Experiment 1. Identified inclusions, no noise, <5 = 0.00. 



Xl 


X2 


X3 


V 


1.640 


-0.510 


0.520 


1.20000 


-1.430 


-0.500 


0.580 


0.50000 


1.220 


0.570 


0.370 


0.70000 


1.410 


0.230 


0.740 


0.61000 


-0.220 


0.470 


0.270 


0.70000 


-1.410 


0.230 


0.174 


0.60000 


Table 3 Experiment 1 


Xl 






V 


1.645 


-0.507 


0.525 


1.24243 


1.215 


0.609 


0.376 


0.67626 


-0.216 


0.465 


0.275 


0.69180 


-1.395 


0.248 


0.177 


0.60747 



. Identified inclusions, S = 0.05. 



Table 4 Experiment 2. Identified inclusions, no noise, <5 = 0.00. 



Xl 


X2 


X3 


V 


1.656 


-0.409 


0.857 


1.75451 


-1.476 


-0.475 


0.620 


0.48823 


1.209 


0.605 


0.382 


0.60886 


-0.225 


0.469 


0.266 


0.69805 


-1.406 


0.228 


0.159 


0.59372 


Table 5 Experiment 2 


Xl 


X2 


£3 


V 


1.575 


-0.523 


0.735 


1.40827 


-1.628 


-0.447 


0.229 


1.46256 


1.197 


0.785 


0.578 


0.53266 


-0.221 


0.460 


0.231 


0.67803 



. Identified inclusions, S = 0.05. 



more difficult identification problem. The program identified less than 6 original 
inhomogeneities in some runs, thus Tables 4 and 5 contain less than 6 lines. As it 
can be seen, only the strongest inhomogeneities have been identified in these cases. 

As in Experiment 1, we have made 10 independent runs of the program for 
each of the noise levels S = 0.00, S = 0.02, and 8 = 0.05. In the noiseless case 
S = 0.00, 5 or 6 original inhomogeneities were identified every time. In 3 runs an 
additional inhomogencity (an artifact) has appeared. For 5 = 0.02 and S = 0.05 
the identification has deteriorated, with 4 inhomogeneities recovered. The number 
of deterministic minimizations and the computational times spent doing them are 
somewhat higher than the ones reported for Experiment 1. 

In general, the execution times were less than 2 minutes on a 333MHz PC. 
As it can be seen from the results, the method achieves a perfect identification in 
the Experiment #1 when no noise is present. The identification deteriorates in the 
presence of noise, as well as if the sources and detectors are not located directly 
above the search area. Still the inclusions with the highest intensity and the closest 
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♦ Sources 

■ Detectors 

O Inclusions 

X Identified Objects 



Figure 2 Inclusions and Identified objects for Experiment 1, 8 = 0.00. X3 
coordinate is not shown. 



ones to the surface are identified, while the deepest and the weakest are lost. This 
can be expected, since their influence on the cost functional is becoming comparable 
with the background noise in the data. 

In summary, the proposed method for the identification of small inclusions can 
be used in geophysics, medicine and technology It can be useful in the development 
of new approaches to ultrasound mammography. It can also be used for localization 
of holes and cracks in metals and other materials, as well as for finding mines from 
surface measurements of acoustic pressure and possibly in other problems of interest 
in various applications. 

The HSD minimization method is a specially designed low-dimensional mini- 
mization method, which is well suited for many inverse type problems. The prob- 
lems do not necessarily have to be within the Born approximation range. It is highly 
desirable to study applications of this method, and to compare its performance to 
other competitive methods. 
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Figure 3 Inclusions and Identified objects for Experiment 2, 8 = 0.00. X3 
coordinate is not shown. 
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